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Abstract 

This  paper  provides  a  theory  for  quantifying  the  hysteresis  and  constitutive  nonlinearities  in¬ 
herent  to  piezoceramic  compounds  through  a  combination  of  free  energy  analysis  and  stochastic 
homogenization  techniques.  In  the  first  step  of  the  model  development,  Helmholtz  and  Gibbs  free 
energy  relations  are  constructed  at  the  lattice  or  domain  level  to  quantify  the  relation  between  the 
field  and  polarization  in  homogeneous,  single  crystal  compounds  which  exhibit  uniform  effective 
fields.  The  effects  of  material  nonhomogeneities,  polycrystallinity,  and  variable  effective  fields  are 
subsequently  incorporated  through  the  assumption  that  certain  physical  parameters,  including  the 
local  coercive  and  effective  fields,  are  randomly  distributed  and  hence  manifestations  of  stochastic 
density  functions  associated  with  the  material.  Stochastic  homogenization  in  this  manner  provides 
low-order  macroscopic  models  with  effective  parameters  that  can  be  correlated  with  physical  prop¬ 
erties  of  the  data.  This  facilitates  the  identification  of  parameters  for  model  construction,  model 
updating  to  accommodate  changing  operating  conditions,  and  control  design  utilizing  model-based 
inverse  compensators.  Attributes  of  the  model,  including  the  guaranteed  closure  of  biased  minor 
loops  in  quasistatic  drive  regimes,  are  illustrated  through  examples. 
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1  Introduction 


The  capability  of  piezoelectric  materials  to  both  actuate  and  sense  derives  from  the  noncentrosym- 
rnetric  nature  of  the  compounds.  At  very  low  input  field  levels,  changes  in  the  ionic  structure 
produce  reversible  changes  in  the  polarization  whereas  dipole  switching  at  higher  field  inputs  yields 
irreversible  increments  in  the  polarization.  This  generates  strains  in  the  material  and  provides  it  with 
actuator  capabilities.  Alternatively,  applied  stresses  also  alter  the  ionic  configuration  which  gener¬ 
ates  the  voltages  measured  in  piezoelectric  sensors.  The  two  mechanisms  are  respectively  termed  the 
converse  and  direct  piezoelectric  effects. 

The  coupled  converse  and  direct  electromechanical  effects  are  highly  sensitive  and  repeatable 
which  makes  PZT  transducers  the  present  choice  for  applications  such  as  nanopositioning  and  sensing. 
For  example,  the  PZT  positioning  elements  in  an  atomic  force  microscope  (AFM)  can  be  used  to 
achieve  angstrom-level  displacements  while  PZT  sensors  are  presently  being  investigated  for  use  in 
multistage  nanopositioners  [7] .  However,  the  polar  mechanisms  which  provide  piezoelectric  materials 
with  the  dual  converse  and  direct  effects,  and  extreme  electromechanical  sensitivity,  also  produce 
varying  degrees  of  hysteresis  and  constitutive  nonlinearities  as  illustrated  in  Figure  1  for  PZT5A. 
Hysteresis  is  an  inherent  property  of  noncentrosymmetric  compounds  at  all  drive  levels  and  is  due  to 
the  irreversible  changes  which  accompany  dipole  switching;  furthermore,  saturation  at  the  domain 
level  and  material  nonhomogeneities  contribute  nonlinear  effects.  Both  the  hysteresis  and  constitutive 
nonlinearities  must  be  accommodated  for  high  performance  applications  utilizing  PZT  actuators  and 
sensors. 

For  a  broad  range  of  applications,  feedback  laws  can  be  employed  to  mitigate  the  deleterious 
effects  of  hysteresis  and  constitutive  nonlinearities.  This  has  led  to  the  successful  use  of  piezoelectric 
transducers  in  applications  ranging  from  hybrid  motor  design  to  structural  acoustic  control  (e.g.,  see 
[4,  9,  14,  38]).  In  other  regimes,  however,  noise  to  signal  ratios  and  fundamental  control  issues  limit 
the  degree  to  which  feedback  design  can  solely  be  employed  to  linearize  the  response  of  piezoelectric 
actuators.  For  example,  the  positioning  elements  in  atomic  force  microscopes  and  nanopositioners 
are  comprised  of  stacked  or  cylindrical  PZT  actuators.  At  low  drive  frequencies,  high  gain  feedback 
laws  are  presently  employed  to  attenuate  hysteresis  and  nonlinearities  thus  leading  to  the  phenomenal 
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Figure  1.  Hysteresis  measured  at  200  mHz  in  PZT5A  for  peak  inputs  of  600  V,  800  V,  1000  V  and 
1600  V. 
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success  of  the  devices.  However,  at  the  higher  frequencies  required  for  large  scale  product  diagnostics 
or  real-time  monitoring  of  biological  processes,  the  efficacy  of  feedback  laws  is  diminished  by  inherent 
thermal  and  measurement  noise.  Robust  control  design  can  be  used  to  extend  the  frequency  ranges 
of  operation  [27],  but  the  loop  shaping  and  gains  required  to  attenuate  hysteresis  have  the  negative 
effect  of  accentuating  high  frequency  noise.  One  technique  for  circumventing  these  limitations  is 
to  develop  feedforward  or  feedback  loops  which  utilize  highly  accurate  and  efficient  model-based 
inverse  compensators.  Models  designed  for  this  use  must  satisfy  at  least  three  properties:  (i)  they 
must  accommodate  transient  dynamics,  (ii)  they  must  guarantee  closure  of  biased  minor  loops,  and 
(iii)  they  must  be  sufficiently  efficient  to  permit  real-time  control  implementation.  In  this  paper, 
we  develop  a  macroscopic  polarization  model  through  the  combination  of  free  energy  principles  at 
the  lattice  level  and  stochastic  homogenization  techniques  which  guarantees  properties  (i)  and  (ii). 
Furthermore,  initial  investigations  attest  to  its  potential  for  real-time  control  implementation  (iii). 

To  provide  a  context  for  the  approach,  we  first  summarize  techniques  that  have  recently  been 
developed  for  quantifying  the  hysteretic  field-polarization  relation  in  piezoelectric  materials.  These 
techniques  can  be  roughly  categorized  as  employing  energy  principles,  phenomenological  principles, 
or  a  combination  of  the  two,  to  construct  microscopic,  mesoscopic  (domain  or  lattice  level),  or 
macroscopic  models. 

Microscopic  and  mesoscopic  models  typically  employ  fundamental  electromagnetic  or  elastic  en¬ 
ergy  relations  to  quantify  the  nonlinear  dependence  of  the  polarization  P  on  input  fields  E  [19]. 
These  theories  provide  a  framework  for  fundamentally  quantifying  material  properties  at  the  lat¬ 
tice  level  and  may  be  necessary  for  optimal  material  design.  However,  for  control  applications,  the 
large  number  of  required  parameters  and  states  precludes  real-time  implementation  using  present 
hardware. 

Macroscopic  models  are  typically  based  on  phenomenological  principles,  thermodynamic  tenets, 
or  energy  formulations  employed  in  concert  with  homogenization  techniques.  The  former  category 
includes  Preisach  models,  which  were  originated  for  magnetic  hysteresis  [22],  and  have  subsequently 
been  extended  to  piezoceramic  materials  [10,  26].  The  advantage  of  Preisach  theory  lies  in  its 
generality  and  strong  mathematical  foundations  which  provide  a  framework  for  quantifying  hysteresis 
when  the  underlying  physics  is  poorly  understood.  However,  the  generality  of  the  technique  also 
yields  models  which  have  a  large  number  of  nonphysical  parameters.  Hence  physical  attributes  of  the 
data  are  difficult  to  utilize  when  identifying  parameters  or  updating  models  to  accommodate  changing 
operating  conditions  although  a  number  of  extensions  to  the  classical  Preisach  theory  have  recently 
been  proposed  to  facilitate  identification  of  parameters  through  correlation  with  physical  principles 
[8,  20].  Moreover,  the  original  Preisach  theory  does  not  accommodate  reversible  effects  or  variable 
temperature,  broadband  operating  conditions,  and  the  modifications  required  to  accommodate  these 
effects  can  significantly  diminish  the  efficiency  of  resulting  models. 

Macroscopic  models  based  on  energy  techniques  provide  a  compromise  between  microscopic  or 
mesoscopic  models  and  solely  phenomenological  models.  Models  in  this  category  include  the  theory 
of  Chen  and  Lynch  [5] ,  quasistatic  hysteresis  models  of  Huang  and  Tiersten  [12]  and  the  domain  wall 
theory  of  Smith,  Horn  and  Ounaies  [32,  33].  While  the  underlying  assumptions,  specific  formulations, 
and  final  goals  differ  in  the  respective  approaches,  similar  strategies  are  employed  when  constructing 
models.  In  each  case,  energy  relations  are  derived  at  the  lattice  or  domain  level  and  averaging 
techniques  are  invoked  to  obtain  macroscopic  models  having  a  small  number  of  effective  parameters. 
For  example,  in  the  theory  of  [5],  energy  relations  at  the  grain  level  are  combined  with  macroscopic 
averaging  over  grain  configurations  to  quantify  strains  and  polarization  in  the  aggregate  material. 
In  the  domain  wall  theory  of  [32,  33],  energy  principles  are  employed  to  quantify  local  changes  in 
polarization  due  to  domain  wall  movement.  These  effects  are  then  averaged  to  obtain  macroscopic 
models  for  bulk  material  characterization  and  control  design. 
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The  theory  presented  here  is  based  on  energy  relations  derived  at  the  lattice  or  domain  level  with 
stochastic  homogenization  techniques  employed  to  construct  macroscopic  models  having  a  small 
number  of  effective  parameters.  In  the  first  step  of  the  development,  Boltzmann  principles  are  used 
to  construct  an  expression  for  the  Helmholtz  energy  through  a  balance  of  the  internal  energy  due 
to  positive  or  negative  dipole  configurations  and  entropy  effects.  The  inclusion  of  the  electrostatic 
work  term  provides  a  Gibbs  relation  which  quantifies  changes  in  the  energy  landscape  due  to  an 
applied  field.  It  is  illustrated  that  minimization  of  the  Gibbs  energy  yields  a  local  polarization  model 
that  quantifies  the  phase  transition  from  the  ferroelectric  to  paraelectric  state  as  temperatures  are 
increased  through  the  Curie  point.  Furthermore,  this  local  relation  is  also  the  Ising  model  employed 
as  an  anhysteretic  kernel  in  the  domain  wall  theory  of  [32,  33].  For  implementation  purposes,  two 
asymptotic  approximations  are  invoked:  (i)  a  quadratic  approximation  to  the  Gibbs  energy  is  con¬ 
structed  for  fixed  temperature  regimes,  and  (ii)  an  algebraic  model  is  constructed  for  the  limiting 
case  of  negligible  thermal  activation.  This  provides  a  highly  efficient  technique  for  quantifying  the 
E-P  relation  in  single  crystal,  homogeneous  materials  with  uniform  effective  fields.  In  the  final  step 
of  the  model  development,  the  effects  of  polycrystallinity,  material  nonhomogeneities,  and  nonuni- 
forrn  effective  fields  are  incorporated  by  considering  physical  parameters  such  as  the  coercive  and 
effective  fields  to  be  manifestations  of  stochastic  distributions.  This  is  motivated  by  the  assumption 
that  different  domains  have  different  energy  characteristics,  and  it  yields  macroscopic  models  with 
parameters  that  effectively  homogenize  or  average  the  material  properties.  It  is  illustrated  through 
fits  to  experimental  data  that  in  spite  of  the  incorporation  of  stochastic  averages,  several  of  the 
effective  parameters  can  be  directly  correlated  with  physical  properties  of  the  data  to  aid  parameter 
identification  and  model  updating.  It  is  also  illustrated  that  the  model  enforces  the  deletion  property 
and  guarantees  closure  of  both  symmetric  and  asymmetric,  biased  minor  loops.  The  model  does  not 
enforce  congruency  near  saturation  which  reflects  the  observation  that  the  measured  E-P  response 
of  the  materials  also  does  not  exhibit  congruency  in  these  regions. 

We  note  that  while  the  specific  motivation  differs,  analogous  concepts  involving  stochastically  dis¬ 
tributed  parameters  are  employed  in  Preisach  formulations  for  magnetic  compounds  [8]  and  discrete 
models  for  shape  memory  alloys  (SMA)  based  on  elastic  chains  constructed  from  bi-stable  elements 
[23,  24,  25]. 

To  place  the  theory  in  a  broader  context,  we  note  that  the  free  energy  analysis  used  to  construct 
the  Helmholtz  and  Gibbs  relations  is  an  extension  of  the  Muller- Achenbach-Seelecke  theory  for  SMA 
[21,  29,  30]  whereas  an  analogous  technique  utilizing  free  energy  relations  in  concert  with  stochastic 
distributions  for  the  coercive  and  effective  fields  has  been  developed  and  implemented  for  ferromag¬ 
netic  compounds  [31].  Hence  the  technique  provides  a  general  methodology  for  quantifying  hysteresis 
and  constitutive  nonlinearities  inherent  to  a  broad  range  of  ferroic  compounds  [35] .  Furthermore,  it 
is  illustrated  in  [37]  that  the  theory  provides  an  energy  basis  for  Preisach  models  with  two  important 
differences:  (i)  the  physical  nature  of  parameters  in  the  proposed  model  facilitates  correlation  with 
properties  of  the  data,  and  (ii)  temperature  and  relaxation  dependencies  are  incorporated  in  the  basis 
rather  than  in  parameters  as  is  the  case  for  Preisach  formulations.  The  latter  property  eliminates  the 
necessity  of  vector-valued  weights  or  lookup  tables  for  material  characterization  and  control  design. 


2  Free  Energies  for  Materials  with  Homogeneous  Lattices 

Energy  formulations  for  commonly  employed  ferroelectric  materials  can  be  motivated  by  changes 
which  occur  in  the  ionic  structure  during  phase  transitions  and  in  response  to  input  fields  and 
stresses.  To  simplify  the  discussion,  we  will  focus  on  barium  titanate  BaTiC>3  and  the  piezoelectric 
compound  Pb(Zr,Ti)C>3  or  PZT.  As  detailed  in  [18],  currently  employed  PZT  transducer  materials 
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are  comprised  of  PbTii_r03  and  PbZrx03  with  x  chosen  to  optimize  electromechanical  coupling. 
The  motivating  discussion  will  emphasize  PbTiC>3  but  analogous  conclusions  hold  for  PbZrC>3  and 
BaTi03. 

These  compounds  are  isostructural  with  the  mineral  perovskite  (CaTiOs)  and  exhibit  what  is 
termed  a  perovskite  structure  comprised  of  a  cubic  form  at  temperatures  T  above  the  Curie  point  Tc 
and  a  tetragonal  form  for  T  <TC.  We  initially  consider  the  idealized  case  of  homogeneous  materials 
having  uniform  lattices;  hence  for  T  >  Tc,  a  unit  cell  at  any  point  in  the  material  will  have  the  cubic 
ionic  structure  illustrated  for  PbTii-^Os  in  Figure  2a.  At  temperatures  below  Tc,  the  materials 
distorts  from  the  cubic  to  tetragonal  form  through  the  biasing  of  Ti4+  ions  toward  02~  pairs  as 
illustrated  for  BaTiC>3  on  page  71  of  [16].  In  the  absence  of  an  applied  electric  field  E,  this  ionic 
configuration  produces  a  double  well  potential  energy  profile  which  varies  as  a  function  of  the  Ti1+ 
position.  As  depicted  in  Figure  3,  the  application  of  an  electric  field  distorts  the  energy  landscape  and 
a  dipole  switch  occurs  when  the  equilibrium  value  determining  the  Ti4+  position  exceeds  the  unstable 
equilibrium  due  to  the  central  02~  pairs.  At  the  macroscopic  scale,  this  produces  a  discontinuous 
jump  in  the  polarization  as  experimentally  illustrated  for  single  crystal  BaTiC>3  on  pages  72-76  of  [16]. 

2.1  Helmholtz  Energy 

We  consider  two  techniques  to  quantify  the  Helmholtz  free  energy  depicted  in  Figure  2c  and  3b; 
the  first  is  based  on  statistical  mechanics  principles  and  the  second,  while  motivated  by  the  first,  is 
phenomenological . 

Based  on  the  assumption  of  material  homogeneity,  we  consider  a  uniform  lattice  of  volume  V 


Figure  2.  (a)  Perovskite  structure  of  PbTiC>3  in  the  cubic  form  above  Tc.  (b)  Tetragonal  structure 
of  PbTi03  for  T  <TC  and  resulting  spontaneous  polarization,  (c)  Helmholtz  free  energy  as  a  function 
of  Ti  position  along  the  ^3-axis. 
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(b) 


Figure  3.  (a)  Helmholtz  energy  ij’  and  Gibbs  energy  G  for  increasing  field  E.  (b)  Local  polarization 
P  as  a  function  of  E  for  a  homogeneous,  isotropic  material. 


and  mass  m  having  N  cells  of  the  form  depicted  in  Figure  2.  Each  cell  is  assumed  to  have  dipole 
orientation  Si  =  ±1  and  dipole  moment  po  so  the  polarization  for  the  lattice  is 

N 

p  = 

(i) 

=  ^(N+-N_). 

Here  Ps  =  Npo/V  denotes  the  saturation  polarization  which  occurs  when  all  dipoles  are  positively 
aligned  and  N+  and  N-  respectively  denote  the  number  of  cells  having  positive  and  negative  orien¬ 
tations.  From  the  second  equality  in  (1)  and  the  fact  that  N+  +  iV_  =  N,  it  follows  that 


JV_  =  —  1  - 


N+  = 


To  compute  the  internal  energy  due  to  dipole  reorientation,  we  let  4>o  denote  the  energy  required 
to  reorient  a  single  dipole  when  the  lattice  is  completely  ordered  ( P  =  Ps).  For  a  general  lattice  or¬ 
dering,  the  energies  required  to  convert  a  dipole  with  positive  orientation  to  negative,  and  conversely, 
are  respectively 

N+  N- 

$+-  =  -J^$0  ,  $-+  =  0 •  (3) 

We  point  out  that  these  energy  expressions  are  derived  under  the  assumption  that  dipoles  interact 
only  with  adjacent  neighbors  (e.g.,  see  [11]). 

The  change  in  the  internal  energy  due  to  dipole  reorientation  can  then  be  expressed  as 


dU  =  [<S>+-dN-  +  $-+dN+}-  .  (4) 

By  utilizing  the  relations  (2)  and  (3),  the  expression  (4)  can  be  integrated  to  obtain  the  internal 
energy  relation 
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where  Uq  denotes  the  energy  for  the  completely  ordered  state.  Since  internal  energy  measures  are 
relative,  we  take  Uq  =  0  in  subsequent  relations. 

The  Helmholtz  energy  for  the  lattice  is  given  by 

i/j  =  U  —  ST 


where  S  denotes  the  entropy.  From  classical  statistical  mechanics  arguments  in  combination  with 
Stirling’s  formula  (e.g.,  see  [11]),  the  entropy  is  given  by 


5  = 


k  i 

v ln 


N 

N, 


k, 
v ln 


N\ 


NJ.  N+\ 


kN 

~V~ 


ln2-i±^lnfl  +  ^ 
2  V  P, 


2  {  R 


-kN 
2  VP, 


P  _L  p 

P  ln  I  )  +  Ps  In  |  1 


P 


+  S0 


where  k  denotes  Boltzmann’s  constant  and  Sq  =  ^ln2.  As  with  Uq,  we  neglect  So  in  the  final 
relation  for  the  Helmholtz  energy  since  we  are  interested  in  a  relative,  rather  than  absolute,  measure 
of  energy. 

The  resulting  Helmholtz  free  energy  is 


*KP,T)  =  %^[1  -(P/P,f]+TkN 


2V 


2  HP, 


AA  [i  -  (p/p,)2]  +  EhT 


2  Tr 


Pin 


Pin 


P  +  Ps 
P,-P 


P  +  Ps 
Ps-P 


+  Ps  ln(l  —  (P/P,)2) 


+  Psln(l  -  (P/Ps)2) 


(5) 


In  the  second  expression  of  (5),  Eh  =  yjr  is  a  bias  field  and  Tc  =  ^  denotes  the  Curie  tempera¬ 
ture  for  the  material.  As  illustrated  in  Figure  4,  the  relation  (5)  yields  a  double  well  potential  at 
temperatures  T  <TC  whereas  behavior  indicative  of  paraelectric  materials  is  reflected  by  the  single 
potential  produced  at  T  >  Tc.  However,  the  logarithmic  nature  of  the  entropic  term  reduces  the 
efficiency  of  algorithms  which  employ  this  relation  and  makes  it  difficult  to  correlate  parameters  in 
the  model  with  physically  measured  properties  of  the  data. 

A  second  technique  for  constructing  the  free  energy,  which  addresses  these  difficulties,  is  based 
on  the  observation  that  for  fixed  temperatures,  a  first-order  approximation  to  the  Helmholtz  relation 
(5)  exhibits  a  quadratic  dependence  on  the  polarization  in  neighborhoods  of  all  three  equilibria.  This 
motivates  the  piecewise  quadratic  definition 

f  ^(P  +  Pr)2 


P(P)  = 


HP-Pr? 

HP!  -  Pr)  (g  -  Pr 


P  <  -  Pi 
P>  Pi 
|P|  <  P/ 


(6) 


for  the  Helmholtz  energy  for  fixed  temperature  regimes.  As  illustrated  in  Figure  3a,  Pj  and  Pr 
respectively  denote  the  positive  inflection  point  and  polarization  at  which  the  minimum  occurs.  The 
relation  between  these  points  and  local  properties  of  the  hysteresis  kernel  will  be  established  in 
subsequent  discussion.  It  will  also  be  illustrated  later  that  rj  can  be  interpreted  as  the  reciprocal 
slope  of  the  hysteron  in  the  post-switching  linear  regime. 
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Figure  4.  Helmholtz  free  energy  specified  by  (5)  for  (a)  T  <  Tc,  and  (b)  T  >  Tc. 

2.2  Gibbs  Energy 

To  construct  a  Gibbs  free  energy  which  exhibits  the  behavior  depicted  in  Figure  3a,  the  relation 
Ue  =  —  p  •  E  quantifying  the  potential  energy  of  a  dipole  p  in  the  field  E  is  combined  with  the 
Helmholtz  energy  throughout  the  lattice  to  yield 

G  =  ip  —  EP  (7) 

for  ip  given  by  (5)  or  (6).  In  the  absence  of  applied  stresses,  the  Gibbs  relation  (7)  is  used  to 
characterize  the  energy  landscape  for  homogeneous  materials  having  a  uniform  lattice  structure. 

2.3  Local  Average  Polarization 

The  probability  of  obtaining  the  energy  level  G  for  a  lattice  volume  V  is  specified  through 
Boltzmann  principles  to  be 

H(G)  =  Ce~GV/kT  (8) 

where  C  is  an  integration  constant  chosen  to  yield  a  probability  of  1  for  integration  over  all  admissible 
input  values.  For  model  identification,  the  volume  V  is  typically  chosen  to  yield  relaxation  behavior 
appropriate  for  the  material  under  consideration. 

For  a  uniform  input  field  E,  the  local  average  polarization  at  the  lattice  level  is  given  by 

P  =  x+  (P+)  +  x-  (P-)  (9) 

where  x+  and  denote  the  fractions  of  dipoles  having  positive  and  negative  orientations  and  (P+) 
and  (P~)  are  the  expected  polarization  levels  associated  with  positively  and  negatively  oriented 
dipoles. 

The  evolution  of  the  dipole  fractions  is  specified  by  the  differential  equations 

x+  =  -p+-x+  +  p-+x - 
X-  =  —p _ \-X-  +  p. | _ x+ 

which  can  be  simplified  to 

x+  =  -P+-X+  +  P-+(1  -  x+)  (10) 
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through  the  identity  x+-\-x~  =  1.  Here  p. ( _ denotes  the  likelihood  of  a  dipole  switching  from  positive 

to  negative  orientation  whereas  p _ |_  denotes  the  likelihood  of  switching  from  negative  to  positive 

(we  avoid  defining  p. \ and  p as  probabilities  since  they  can  be  unbounded).  The  likelihoods  are 

computed  by  specifying  the  probability  of  achieving  the  energy  required  for  a  jump  multiplied  by 
the  frequency  at  which  jumps  are  attempted.  This  yields  the  relations 


P+-  = 


P-+  = 


kT 


o-G(E,P0(T),T)V/kT 


e-G(E,P,T)V/kTdp 


2irmV2/3  f°° 

JPo(T ) 

-G(E,-P0(T),T)V/kT 


kT 


2irmV2/3  rp b(T) 


-G(E,P,T)V/kT  dp 


where  Po(T)  is  the  unstable  equilibrium  of  G  and  m  denotes  the  mass  of  the  lattice.  The  denominator 
in  both  cases  arises  from  evaluation  of  the  integration  constant  C.  When  implementing  the  model, 
we  typically  replace  Pq(T)  by  the  inflections  points  Pj  and  —Pi,  respectively,  in  the  relations  for 
p. ( _ and  p _ to  obtain 


P+-  = 


kT 


0—G(E,PI,T)V/kT 


2irmV2/3 


-G(E,P,T)V/kT  dp 


'Pi 


P-+  = 


kT 


a—G(E,—PI,T)V/kT 


2irmV2/3 


'—Pi 


e-G(E,P,T)V/kTdp 


(11) 


This  simplifies  the  approximation  of  the  integrals  and  can  be  motivated  by  observing  that  if  one 
considers  the  forces  due  to  the  applied  field,  maximum  restoring  forces  occur  at  Pj  and  —Pi 
(e.g.,  see  pages  332-333  of  [6]).  Furthermore,  for  materials  with  low  thermal  activation,  the  points 
Po  and  —Pi  coincide  in  the  limit  as  thermal  activation  is  reduced  to  zero  for  positive  fields  while  Pj 
and  Po  coincide  for  negative  fields  as  illustrated  in  Figure  5. 

The  expected  polarizations  are  given  by 

roo  fpo(T) 

(P+)=  Pp(G)dP  ,  <P_>=/  Pp(G)dP. 

JP0{T)  J- oo 


Specification  of  the  probabilities  using  (8)  and  formulation  of  the  integrals  in  terms  of  the  inflection 
points  yields  the  relations 


/  Pe—G(E,P,T)V/kT  dp 

(P+)  =  ~^So -  ,  (P-)  = 

/  e-G(E,P,T)V/kTdp 
JPi 


—  P, 


Pe-G(E,P)V/kTdp 


r-Pi 


-G(E,P)V/kTdp 


(12) 


quantifying  the  expected  polarizations  respectively  due  to  positively  and  negatively  oriented  dipoles. 

The  summed  products  of  the  expectations  and  phase  fractions  (9)  quantifies  the  local  polarization 
P  at  the  lattice  level.  For  materials  which  are  homogeneous  and  isotropic,  this  local  polarization 
will  be  uniform  throughout  the  material  and  hence  it  will  also  quantify  the  macroscopic  polarization 


(a)  (b) 


Figure  5.  (a)  Gibbs  energy  profile  with  high  levels  ( - )  and  low  levels  ( - )  of  thermal  activation 

in  the  Boltzmann  distribution  /x(G)  =  Ce~GV^kT .  (b)  Local  polarization  P  given  by  equation  (9) 

with  high  thermal  activation  ( - )  and  limiting  polarization  P  specified  by  (22)  in  the  absence  of 

thermal  activation  ( - ). 


generated  in  response  to  an  input  field.  With  reasonable  accuracy,  this  model  could  be  used  to 
quantify  the  single  crystal  BaTiC>3  behavior  depicted  on  pages  72-76  of  [16].  Extensions  to  the 
model  to  accommodate  nonhonrogeneous  material  and  field  attributes  will  be  provided  in  Section  3. 

The  probabilistic  nature  of  dipole  switching  produces  the  gradual  transitions  depicted  in  Fig¬ 
ures  3b  and  5b,  with  the  degree  to  which  transitions  are  mollified  being  dependent  on  the  ratio 
between  GV  and  kT  in  the  Boltzmann  relation  (8).  Large  values  of  kT,  relative  to  GV,  characterize 
regimes  in  which  thermal  activation  is  prominent  which  in  turn  produces  smoother  transitions  since, 
for  a  fixed  field  level,  dipoles  have  a  higher  probability  of  achieving  the  thermal  energy  required  to 
overcome  energy  barriers. 

For  materials  and  operating  regimes  in  which  the  relaxation  effects  due  to  thermal  activation  are 
negligible,  the  differential  equation  model  which  yields  (9)  can  be  simplified  to  a  purely  algebraic 
model  through  asymptotic  analysis.  This  can  improve  the  efficiency  of  model-based  characterization 
and  control  algorithms  and  highlight  additional  properties  of  the  local  model. 

As  a  prelude  to  this  asymptotic  analysis,  we  establish  the  following  theorem. 

Theorem  1.  Let  f  be  continuous  on  the  interval  [— L,L\  and  let  {<fj}  be  a  sequence  satisfying  the 
following  properties: 


(i)  4>j  >  0  for  all  j 

(ii)  f  (fj(y)dy  =  1  for  all  j 

(in)  Given  e,5  >  0,  there  exists  jo  such  that 


f>j{y)\f{y)\dy  <  e/3 


for  all  j  >  j0 ■ 
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Then  for  x  6  [— L,  L\,  <pk  *  f  converges  to  f;  that  is 


f  <t>j{x-y)f(y)dy  f(x). 


Proof.  From  (ii)  it  follows  that 

f(x)  =  f(x)  j  <t>j(y)dy  =  j  f(x)</>j(y)dy 

so  that  for  x  G  [—L,L], 

<t>j*f{?)-f{x)  =  J  <j)j{y)f(x-y)dy-  J  <j>j(y)f(x)dy 

=  J  4>j(y)[f(x  -y)~  f{x)}dy . 

Furthermore,  from  the  continuity  of  /,  it  follows  that  for  fixed  e,  there  exists  5  such  that 

I  f(x-y)  -  f(x)  |  <  e/3 

for  |  y  |  <  5.  For  this  5,  we  write 

*  f(x)  ~  fix)  =  [  f>j{y)[f{x  -y)-  f(x)\dy  +  [  </>j{y)[f{x  -  y)  -  f{x)}dy  <  e 

J\y\<5  J\y\>6 

for  sufficiently  large  j.  The  integral  over  the  region  |y|  <  5  is  bounded  by  s/3  due  to  the  continuity 
of  /  whereas  the  integral  over  |y|  >  5  is  bounded  by  2e/3  for  sufficiently  large  j  as  a  result  of  (iii). 
The  convergence  follows  directly.  □ 

We  note  that  if  we  replace  property  (iii)  by  the  requirement 


/  <~h(y)dy  <  £ 

d\y\>S 

and  add  the  assumption  that  /  is  measurable  and  bounded  on  ]R,  the  sequence  (j>j  is  termed  a  Dirac 
sequence  on  ]R,  and  Theorem  1  is  a  1-D  version  of  Theorem  3.1  from  page  228  of  [13]. 


2.4  Asymptotic  Polarization  Relation  in  Absence  of  Thermal  Activation 

We  now  consider  initial  dipole  fractions  x~,  x+,  and  a  positive  field  E  for  which  G{Pmin)  <  G{Po ) 
and  G{Pmin)  <  G(Pmin )  as  depicted  in  Figure  5a.  For  simplicity,  we  consider  the  piecewise  quadratic 
Helmholtz  free  energy  model  (6)  and  note  that  analogous  asymptotic  analysis  holds  for  the  statistical 
mechanics  model  (5).  In  this  case,  the  minima 


Pmin{E )  = - Pr 

V 


Pmin{E )  = - 1-  Pr 

V 


(13) 


result  from  the  necessary  condition 

by 


(]p  =  (].  The  coercive  field  for  which  Pmin  =  Pq  =  —  Pi  is  given 


Ec  =  v(Pr  ~  Pi)  • 


(14) 
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We  illustrate  first  that  for  kT /V  >  0,  the  Boltzmann  probability  (8)  exhibits  Gaussian  behavior 
in  neighborhoods  of  Pmin  and  Pm,in  with  variance  proportional  to  kT /V.  We  will  also  illustrate  that 
in  the  limit  kT /V  — >  0,  l-i(G)  converges  to  a  Dirac  distribution. 

For  P  <  —Pi,  h{G)  can  be  formulated  as 


MG) 


a-[h l{P+PR)2-EP\V/kT 


r -Pr 


-[\n{P+PRf-EP]V/kT  dp 


e-(P-Pmin)2vV/2kT 


'-P, 


a—(P—Pmin)2,nV/2kT  dp 


C(T,/3)e  (P  Prnin)2/2P2 


(15) 


where 


Hw 


C(T,(3)  = 


(16) 


r-Pr 


,-{P-Pmin)2/202(T)dp 


-1 


We  now  let  j  =  1//3  and  define  the  sequence 


f  C(T,j)e-(p-p-^'2/2  ,  P<-Pj 

-  \  o  ,  p  >  -Pj . 


Since  satisfies  (i)-(iii)  of  Theorem  1,  and  hence  constitutes  a  Dirac  sequence,  it  follows  that 


lim  u(G)  =  lim  d>j(P) 

kT/V— 0  oo 

=  S(P  —  Pmir 


(17) 


Analogous  behavior  is  exhibited  at  Pmin- 

The  Gaussian  behavior  of  /j,  quantified  by  (15),  is  depicted  in  Figure  5a.  From  the  definition  of 

(3 ,  it  follows  that  a  decrease  in  thermal  activation  is  reflected  as  a  decrease  in  the  variance.  This 

implies  that  a  smaller  number  of  dipoles  achieve  the  energy  required  to  overcome  the  energy  barrier 
which  yields  steeper  transitions  in  the  relation  between  E  and  P  as  depicted  in  Figure  5b. 

We  now  illustrate  that  the  solution  to  the  model  (9)  converges  to  the  piecewise  linear  kernel 
specified  by  (13).  It  is  first  noted  that  as  kT /V  — ►  0,  the  limiting  solution  to  (10)  is 

x+  =  x+  +  (1  —  x+)h(P  +  Pi) 

I  x+  ,  E<EC  (18) 

\  1  ,  E>  Ec 

where  the  local  coercive  field  Ec  is  given  by  (14)  and  h  denotes  the  Heaviside  function.  Corresponding 
values  for  x_  are  determined  through  the  relation  x+  +  =  1. 


11 


The  expected  polarization  due  to  positively  oriented  dipoles  is 


(P+) 


lim 

kT/V^O 


/» OO 

/  Pe~G^v/kTdP 
Jpi 


-G(E,P)V/kT  dp 


’pi 


lim 

j^oo  J 


/*°o 

— /  Pe-(p-prnin)2j2/ 2  Jp 

\/2 pJPi 


\/27T 


•  poo 

! /  e-{P-Pmin)2j2/Z  dp 

27 r  Jpr 


Prr 


(19) 


with  the  limits  in  the  numerator  and  denominator  evaluated  using  Theorem  1  with  /(P)  =  P  and 
/(P)  =  1  for  P  >  Po,  /(P)  =  0  for  P  <  Po,  respectively.  Similarly,  the  limiting  value  of  (P_)  is 

<P_)  =  Pmin.  (20) 

From  (18)-(22),  it  follows  immediately  that  for  the  initial  dipole  distribution  x+  and  X-,  the  local 
polarization  is  given  by 


P{E) 


X—Pmin(P)  T  %+Pmin(E)  0  <  E  <  Ec 

Pmin{E )  E  >  Ec 


(21) 


for  positive  fields  E  as  depicted  in  Figures  3  and  5.  Analogous  results  hold  for  negative  fields. 

To  accommodate  multiple  transitions,  the  local  polarization  resulting  from  the  Helmholtz  relation 
(6)  with  no  thermal  activation  can  be  formulated  using  Preisach  notation  (e.g.,  see  [2,  37])  as 


[p(E-Ec,om  =  i 


[p(E-Ec,om 
e~Pr 
#  +  PR 


where 


[P(P;Pc,<e)](0)  =  ^ 


,  r(f)  =  0 

,  r(t)  /  0  and  P(maxr(t))  =  —  Ec 
,  r(t)  /  0  and  F(maxr(t))  =  Ec 


(22) 


'  f -Pr 


f  +  Pr 


,  E(0)  <  Ec 
,  -Ec  <  E( 0)  <  Ec 
,  E( 0)  >  Pc 


denotes  the  dipole  orientation  yielding  the  initial  polarization.  A  depiction  of  representative  initial 
polarization  values  is  provided  in  Figure  8.  The  transition  times  are  designated  as 


r(t)  =  {t  e  (0,  Tf  \  |  E(t)  =  -Ec  or  E(t)  =  Ec} 

where  Tf  denotes  the  final  time  under  consideration.  The  dependence  of  P  on  the  local  coercive  field 
Ec  defined  in  (14)  is  explicitly  indicated  since  the  parameter  is  considered  to  be  distributed  in  the 
next  section.  We  also  note  that  the  piecewise  linear  models  (21)  or  (22)  can  be  obtained  directly  from 
the  necessary  condition  fp  =  0;  the  asymptotic  analysis  illustrates  the  consistency  of  this  condition 
in  the  limiting  behavior  of  the  full  model  (9). 
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Figure  6.  (a)  Kernel  provided  by  (22),  and  (b)  kernel  provided  by  (23). 


Similar  analysis  can  be  used  to  specify  the  hysteresis  kernel  which  results  when  G  is  constructed 
using  the  Helmholtz  relation  (5)  derived  through  statistical  mechanics  arguments.  The  determination 
of  P  through  solution  of  =  0  in  this  case  yields 


P(E-  Eh) 


Pstanh 


(E  +  EhP/Ps\ 
V  EhT/Tc  ) 


Pstanh 


E  +  aP\ 
a(T)  ) 


(23) 


where 


(24) 


The  behavior  of  P  specified  by  (23)  is  compared  in  Figure  6  with  that  defined  by  (21)  or  (22)  which 
was  derived  from  the  piecewise  quadratic  Helmholtz  energy.  While  the  two  representations  yield 
similar  kernel  behavior  as  dipole  switching  occurs,  the  model  (23)  predicts  a  local  saturation  value 
of  Ps  and  decreasing  slope  for  increasing  field  whereas  the  model  (22)  predicts  linear  behavior  after 
dipole  switching  with  slope 

The  Ising  model  (23)  has  been  employed  in  a  number  of  hysteresis  models  for  ferroelectric  mate¬ 
rials  and  its  ubiquity  is  due  to  the  common  underlying  assumption  that  dipoles  are  aligned  in  only 
two  configurations:  in  the  direction  of  the  applied  field  or  diametrically  opposite  to  it.  As  detailed  in 
[32],  quantification  of  the  electrostatic  energy  under  this  assumption  yields  the  Ising  model  whereas 
the  relaxed  assumption  that  dipoles  can  orient  uniformly  yields  a  Langevin  model,  which  agrees  with 
the  Ising  model  through  first-order  terms  but  predicts  a  slower  saturation  rate.  In  both  cases,  these 
models  were  employed  to  quantify  the  anhysteretic  kernel  as  an  initial  step  in  the  development  of  a 
macroscopic  hysteresis  model  [32,  33].  The  theory  developed  here  differs  from  the  domain  wall  theory 
in  [32,  33]  in  the  sense  that  the  Ising  model  directly  quantifies  the  energy  landscape  at  the  lattice 
level  whereas  it  provides  only  an  intermediate  step  in  the  domain  wall  theory.  The  saturation  be¬ 
havior  of  the  Ising  relation  has  also  motivated  its  use  in  phenomenological  models.  Translates  of  the 
form  P  =  Pstanh(F  ±  Ec )  were  employed  by  Zhang  and  Rogers  [39],  and  r{x )  =  tanh(.x)  provides  a 
suitable  choice  for  the  ridge  functions  employed  in  generalized  Preisach,  or  Krasnolselskii-Pokrovskii, 
characterizations  [2,  3].  The  model  presented  here  differs  in  that  it  focuses  on  an  energy  formulation 
for  domain  processes.  This  yields  a  low-order  model  which,  as  illustrated  in  the  examples,  ensures 
closure  of  biased  minor  loops. 
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3  Polycrystalline  Materials  with  Variable  Effective  Fields 


The  local  polarization  relations  (9),  (22)  and  (23)  were  derived  under  the  assumption  of  homogeneous 
and  isotropic  material  properties  and  uniform  effective  fields  Ee  =  E  throughout  the  materials.  The 
consideration  of  the  local  relations  throughout  the  material  yields  global  models  which  predict  instan¬ 
taneous  transitions  at  the  coercive  field  as  illustrated  in  Figures  4,  5  and  6.  While  such  global  models 
can  accurately  quantify  single  crystal  behavior  of  the  type  experimentally  measured  for  BaTiC>3  (e.g., 
see  page  76  of  [16]),  they  do  not  accurately  predict  the  more  gradual  transition  through  the  remanent 
polarization  measured  in  polycrystalline  ferroelectric  materials.  In  this  section,  we  incorporate  the 
effects  of  nonunifornr  lattice  configurations,  polycrystallinity,  and  variable  effective  fields  to  provide 
a  macroscopic  model  which  accurately  characterizes  hysteresis  in  a  variety  of  ferroelectric  materials 
and  ensures  closure  of  biased  minor  loops. 

3.1  Distributions  in  Remanence  Polarization,  Lattice,  and  Coercive  Field 

As  illustrated  in  Figure  7,  nonuniformities  in  the  lattice  produce  a  distribution  of  Helmholtz  and 
Gibbs  free  energy  profiles  which  can  be  manifested  as  variations  in  the  local  coercive  field  and  local 
remanent  polarization  and  can  produce  differing  saturation  behavior  after  dipole  switching.  Addi¬ 
tional  variations  in  the  free  energy  relations  will  be  produced  by  stress  nonhomogeneities,  nonunifornr 
lattice  orientations  across  grain  boundaries,  and  crystalline  anisotropies. 


Figure  7.  (a)  Nonunifornr  lattice  and  polycrystalline  structure  for  PZT;  (b)  Free  energies  associated 
with  lattice  widths  (i)  and  (ii);  (c)  Variations  in  hysteresis  kernel  due  to  differing  free  energies. 
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For  the  piecewise  quadratic  Helmholtz  model  (6),  the  variability  in  the  lattice  structure  can 
be  incorporated  by  considering  Pr,  Pj  or  Ec  =  ij(Pr  —  Pj)  to  be  manifestations  of  an  underlying 
distribution  rather  than  fixed  values  as  assumed  in  the  previous  section  for  single  crystals  having 
uniform  lattices.  For  this  model,  we  consider  variations  in  the  local  coercive  field  Ec  and  make  the 
assumption  that  it  can  be  modeled  by  a  log-normal  distribution  to  incorporate  the  requirement  that 
Ec>  0.  This  yields  the  relation 

aoo 

[P(E)}(t)=  /  \P{E-Ec,Z)}{t)f(Ec)dEc 
Jo 

for  the  macroscopic  polarization  where  the  density  /  is  specified  by 


f(Ec)  =  d  exp 


In  {Ec/Er'2' 
2c 


(25) 


and  P  is  given  by  (9)  or  (22).  Here  c,  ci  and  Ec  are  positive  constants.  In  [8],  it  is  illustrated  that 
if  c  is  small  compared  with  Ec,  the  mean  and  variance  have  the  approximate  values 

(Ec)  ~  Ec  ,  a  k,  2EC  c .  (26) 

As  will  be  illustrated  in  the  examples  of  Section  4,  the  relation  (26)  can  be  used  to  obtain  initial 
parameter  estimates  based  on  attributes  of  measured  experimental  data.  For  materials  in  which  the 
transition  during  hysteresis  is  relatively  gradual,  a  second  choice  for  /  is  a  normal  distribution  with 
mean  Ec  and  variance  a.  The  lower  integration  limit  of  0  should  be  retained  to  enforce  nonnegative 
local  coercive  fields. 


3.2  Nonuniform  Effective  Fields 

The  second  extension  employed  to  obtain  a  macroscopic  model  for  the  polarization  entails  the 
consideration  of  effective  fields  in  the  material.  As  noted  in  [1,  15,  32,  33],  the  applied  field  in 
ferroelectric  materials  is  augmented  by  fields  generated  by  neighboring  dipoles  which  produce  non- 
homogeneous  effective  fields  in  the  materials.  This  produces  variations  about  the  applied  field  which 
can  significantly  alter  the  measured  polarization.  To  incorporate  these  field  variations,  we  consider 
the  effective  field  to  be  normally  distributed  about  the  applied  field.  For  fixed  Ec,  the  polarization 
in  this  case  is  given  by 

/OO 

c2[P(Ee-1Ec,0}(t)e-^E-E^/bdEe.  (27) 

-OO 

The  introduction  of  variations  in  the  effective  field  produces  domain  switching  in  advance  of  the 
remanence  point  in  accordance  with  observations  from  experimental  data. 

The  complete  macroscopic  polarization  model  for  nonhomogeneous,  polycrystalline  materials  with 
variable  effective  fields,  as  based  on  the  piecewise  quadratic  Helmholtz  model  (6),  is  then  given  by 

/*oo  /»oo  _ 

[P(E)](t)  =  C  /  [P{Ee  +  E,Ec,OKt)e-E^be-^Ec/E^^2dEedEc  (28) 

where  P  is  specified  by  (9)  or  (22).  We  note  that  while  the  model  (28)  incorporates  certain  relax¬ 
ation  mechanisms,  it  does  not  incorporate  dynamic  elastic  effects.  Hence  this  formulation  of  the 
polarization  model  should  be  restricted  to  low  frequency  drive  regimes. 
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Similar  arguments  can  be  employed  to  construct  a  macroscopic  polarization  model  based  on  the 
Helmholtz  free  energy  relation  (5)  obtained  through  statistical  mechanics  analysis.  In  this  case,  we 
assume  that  the  bias  field  Eyt  is  a  manifestation  of  an  underlying  log-normal  distribution  to  obtain 
the  global  relation 

/*oo  /*oo  _ 

[P(E)\(t)=C  /  [P(Ee  +  E,Eh,0}{t)e-E^be-^Eh/Eh^2c^dEedEh.  (29) 

J  0  j  —  DO 

Here  P  is  specified  by  (9)  or  (23). 

The  polarization  behavior  predicted  by  (29)  differs  from  that  of  (28)  primarily  at  high  input 
field  levels.  The  E-P  behavior  predicted  by  (28)  reflects  the  linear  behavior  associated  with  the 
hysteresis  kernel  (22)  after  completion  of  dipole  switching  whereas  the  E-P  curve  produced  by  (29) 
saturates  to  zero  slope  due  to  the  behavior  of  the  hyperbolic  tangent  kernels.  While  both  models  are 
appropriate  for  a  number  of  materials,  the  saturation  behavior  and  ease  with  which  the  respective 
models  can  be  implemented  are  effective  criteria  for  choosing  between  the  models  (28)  and  (29)  for 
a  given  application. 

3.3  Implementation  Issues 

To  implement  the  models  (28)  or  (29),  it  is  necessary  to  approximate  the  integrals.  This  can 
be  accomplished  on  the  semi-infinite  domain  using  Gauss-Laguerre  quadrature  and  on  the  infinite 
domain  using  Gauss-Hermite  quadrature  (e.g.,  see  pages  698-699  of  [40]).  Alternatively,  the  expo¬ 
nential  decay  of  the  kernels  can  be  employed  to  truncate  the  domains  to  finite  intervals  appropriate 
for  Gauss-Legendre  formulae  (see  Figure  8a).  In  all  cases,  approximation  of  (28)  yields  expressions 
of  the  form 

Ni  Nj  _  2  - 
[P(E)}(t)  =  C^2^2[P(Eej  +E;F;Ci,e0](Oe'^/V[ln^/^)/2c]2^  (30) 

i= 1  j= 1 

with  a  similar  relation  resulting  from  the  approximation  of  (29).  Here  Ee.,ECi  denote  the  abscissas 
associated  with  respective  quadrature  formulae  and  Vi,u)j  are  the  respective  weights. 

For  the  examples  reported  in  Section  4,  we  employed  composite  Gauss-Legendre  quadrature  on 
truncated  intervals  chosen  to  be  commensurate  with  nontrivial  values  of  the  integrands.  To  illustrate, 
we  provide  the  abscissas  and  weights  employed  in  the  approximation  of  the  integral  (27)  on  the 
truncated  domain  [— L,L]  using  a  four  point  composite  quadrature  rule.  Letting  hj  =  —L  +  jh,  the 
quadrature  points  and  weights  on  each  subinterval  [hj- 1,  hj]  are 

1  \/ 15+2V30  49 h 

2  2V35  ’  12(18+v/30) 

1  y/15-2V30  49/i. 

2  2V35  ’  i2  12(18— \/30) 

1  .  \J  15— 2V30  49 h 

2  ~l~  2^35  ’  J’3  12(18— x/30) 

^  1  \/  15+2y/30  49/i 

2  2^35  ’  J'4  _  12(18+v/30)  ' 

The  quadrature  points  and  initial  polarization  values  ^  for  E  =  0,  P  =  0,  with  the  hysteresis  kernel 
(22),  are  depicted  for  Nq  =  2  (Nj  =  8)  in  Figure  8b. 


Ee.  y  1  —  h  j  -  1  T  h 
Eej  2  —  h  j  —  1  h 
Eej3  —  h  j  —  1  -f-  h 
hhe:j  j  —  hjj  —  1  T  h 
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Figure  8.  (a)  Decay  exhibited  by  f-2{Ee)  =  e  E^b  and  truncated  domain  [—L,L\.  (b)  Gaussian 
quadrature  points  •  and  initial  local  polarization  values  £j  (indicated  by  x)  for  Nj  =  8.  (c)  Log¬ 
normal  density  f(Ec )  =  cie_[ln(Ec/Ec)/2cl2  having  mean  Ec.  (d)  Distribution  of  hysteresis  kernels 
having  coercive  fields  Ec. 


A  second  implementation  issue  concerns  the  manner  through  which  the  conditional  definitions 
in  (22)  are  evaluated  if  this  kernel  is  employed  in  the  hysteresis  model.  While  it  is  algorithmically 
straightforward  to  implement  these  conditions  using  the  transition  times  r(t),  this  must  be  done 
for  all  quadrature  points  Eej  and  ECi  for  each  input  field  value.  Implementation  in  this  manner 
significantly  diminishes  the  speed  of  the  model  evaluation  and  would  likely  prohibit  the  use  of  the 
model  for  real-time  control  design  and  implementation. 

An  alternative  is  to  formulate  the  local  polarization  (22)  as 

P  =  -  +  Pr  A 
V 

where  A  =  1  if  evaluating  on  the  upper  branch  of  the  kernel  and  A  =  —  1  if  evaluating  on  the  lower 
branch.  The  crux  of  the  algorithm  focuses  on  the  efficient  construction  of  A  which  accommodates 
the  vector  nature  of  the  effective  field  values  and  coercive  fields  resulting  from  the  quadrature  of 
the  integrals.  Intuitively,  the  local  polarization  values  associated  with  each  effective  field  value  Eej 
will  jump  when  they  reach  a  coercive  field  ECi.  Because  both  the  effective  and  coercive  field  values 
are  distributed,  as  depicted  in  Figure  8,  this  leads  to  N  x  Nj  relations  which  must  be  checked  for 
each  input  field  value  E.  Hence  for  the  distributed  algorithm,  A  is  an  Nz  x  Nj  matrix  in  which 
the  ij th  element  specifies  whether  the  jth  effective  field  value  Ee.  has  crossed  the  ith  coercive  value 
ECi  to  determine  whether  the  associated  polarization  value  is  on  the  upper  or  lower  branch  of  the 
hysteron.  The  high  efficiency  of  the  algorithm  is  achieved  by  employing  algebraic  matrix  operations 
to  construct  A  rather  than  conditional  statements  implemented  through  if-then  constructs. 


(31) 
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To  illustrate  the  algorithm  used  to  compute  P  using  the  approximate  relation  (30),  with  P  given 
by  (22),  we  first  define  the  following  matrices 


'  -l  •• 

•  -l 

i  • 

•  •  i  ' 

1 _ 

1 - 

^ init  — 

-l  •• 

•  -l 

i  •  ■ 

•  •  i 

; 

NiXNj 

ii 

-  E°Ni  ’ 

£  ■" 

j? 

1 _ 

Ek  +  Eei 

•  •  Efc  +  EeN . 

'  1 

••  1  " 

£k  = 

,  o  = 

Ek  +  Ee  i 

•  •  Ek  +  EeN . 

ivj  j 

NiXNj 

1 

••  1 

NiXNj 

and  weight  vectors 


WT  = 


w\e 


-Eljb 


~ElN./b 

,wNje 


lxNi 


VT  = 


Vie-MEcl/Ec)/2c]^  ...  ?  e 


—  [ln(ECJV  /Ec)/2c]2 


J  lx  TV; 


where  Ej-  =  E{tk)  is  the  kth  value  of  the  input  field.  The  polarization  Pp;  ~  P(Ek)  is  then  specified 
by  the  following  algorithm. 


Algorithm  1. 

A  — 

P=\£c  +  PrO 

for  k  =  1  :  Nk 

P  =  \£k  +  Pr A 

dE  =  Ek  —  Ek- 1 

A  =  sgn((£fc  -  sgn (dE)£c).  *  (P  -  sgn (dE)P).  *  P ) 

P  =  x-£k  +  PrA 

Pk  =  CVTPW 

end 

In  this  algorithm,  .*  indicates  componentwise  matrix  multiplication  and  sgn  denotes  the  signum 
function.  Depending  on  the  methods  used  for  programming,  the  use  of  Algorithm  1  rather  than 
utilizing  conditional  if-then  constructs  can  reduce  runtimes  by  factors  in  excess  of  100  so  that  full 
multiloop  model  simulations  run  in  the  order  of  seconds  on  a  workstation.  This  level  of  efficiency  is 
necessary  to  achieve  real-time  implementation  of  control  algorithms  utilizing  this  model. 


4  Model  Validation 

To  illustrate  attributes  of  the  model,  we  consider  two  sets  of  examples.  In  the  first,  the  capability 
of  the  model  to  characterize  symmetric,  major  loop  properties  of  PZT5A,  PZT5H  and  PZT4  is 
illustrated  through  comparison  with  experimental  data.  For  each  compound,  model  parameters 
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are  estimated  through  a  least  squares  fit  to  data  at  high  drive  levels  and  the  resulting  model  is 
used  to  predict  the  model  behavior  in  lower  drive  regimes.  In  all  cases,  quasistatic  drive  regimes 
are  considered.  The  second  example  illustrates  numerically  the  capability  of  the  model  to  quantify 
biased  minor  loop  behavior  including  Raleigh  loops  at  low  drive  levels  and  multiply  nested  loops  at 
intermediate  levels.  In  concert,  these  examples  illustrate  the  flexibility  of  the  model  for  a  variety  of 
materials  and  operating  conditions. 

4.1  Determination  of  Parameters 

The  continuous  model  (28)  or  discretized  model  (30)  contains  the  parameters  PR,rj,b,Ec,c  and 
C  which  must  be  specified  when  quantifying  a  specific  PZT  compound.  Similar  parameters  arise  in 
the  model  derived  through  statistical  mechanics  principles.  Asymptotic  analysis  can  be  employed  to 
obtain  initial  parameter  choice  which  can  then  be  employed  in  various  least  squares  formulations  to 
determine  parameters  which  optimize  model  fits  and  predictions. 

As  illustrated  in  Figure  5,  Pr  denotes  the  local  remanence  polarization  for  a  domain  and  i/  is  the 
reciprocal  of  the  slope  in  the  E-P  relation  after  switching.  The  inclusion  of  polycrystallinity,  variable 
effective  fields  and  material  nonhomogeneities  through  the  density  analysis  in  Section  3  makes  it 
difficult  to  correlate  the  remanence  polarization  measured  for  the  bulk  material  with  the  local  value 
Pr',  hence  Pr  is  typically  estimated  solely  through  a  least  squares  fit  to  the  data.  Moreover,  for 
the  linear  relations  (22)  or  (31)  for  the  kernel  P,  Pr  and  C  produce  analogous  scaling  in  the  bulk 
polarization  so  that  they  can  be  combined  when  estimating  parameters.  The  slope  of  the  local 
kernel  relations  scales  through  the  stochastic  homogenization  process  so  bulk  measurements  of  the 
reciprocal  slope  (jp  provide  initial  estimates  for  ?/.  The  distribution  /,  defined  in  (25),  quantifies 


Figure  9.  (a)  Asymptotic  relations  between  the  parameters  r/,  Ec,  b  and  the  slope  of  the  P-E  rela¬ 
tion  after  switching,  the  coercive  field  and  the  point  where  switching  commences  before  remanence. 

(b)  The  absolute  metric  (32)  and  Euclidean  metric  (34)  between  data  and  the  model  prediction. 

(c)  Cosine  filter  (33)  used  to  minimize  least  squares  differences  in  the  switching  region. 
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the  coercive  properties  of  the  E-P  relation.  The  mean  Ec  quantifies  the  point  at  which  the  primary 
switching  occurs  as  illustrated  in  Figure  9a;  hence  Ec  is  asymptotically  given  by  the  coercive  field 
for  the  bulk  material.  The  variance  a  ~  2 Ecc  quantifies  the  variability  in  local  coercive  fields  so  that 
materials  with  steep  coercive  transitions  yield  small  values  of  c  whereas  large  parameter  values  are 
required  when  characterizing  materials  with  gradual  bulk  switching.  The  parameter  b  quantifies  the 
variance  in  the  effective  field  which  determines,  in  part,  the  degree  to  which  switching  occurs  before 
remanence  is  reached.  Materials  with  nearly  linear  E-P  relations  at  remanence  yield  small  values  of 
b  whereas  large  values  are  required  to  accommodate  significant  switching  before  remanence. 

Hence  physical  properties  of  the  data  yield  initial  estimates  for  ij  and  Ec  and  provide  quali¬ 
tative  techniques  for  ascertaining  b,  c,  Pr  and  C.  This  significantly  facilitates  the  implementation 
of  least  squares  techniques  used  to  determine  model  parameters  and  update  these  parameters  to 
accommodate  slowly  changing  operating  conditions. 

Three  least  squares  techniques  have  been  considered  to  accommodate  the  switching  and  saturation 
behavior  inherent  to  hysteresis  data.  To  formulate  the  least  squares  problems,  let  ( Ei,Pi),i  = 
1,  •  •  • ,  jV,  denote  the  field  and  corresponding  polarization  data  measured  throughout  the  hysteresis 
cycle.  Furthermore,  let  P(Ep,  q)  denote  parameter-dependent  model  solutions  provided  by  (28),  (29) 
or  (30).  For  admissible  parameters  q  €  Q,  we  then  consider  the  following  optimization  problems: 


JX 

min  y  P(E?;;g)  -  j 

q  7=1 

^  2 

P7  , 

(32) 

M 

min  V'  m 

q  i 

1=1 

P(Ei\q)-Pi 

/J,;  =  cos 

0  Ei  Emin 

2vr-  +a, 

J-^max  Eli  min 

(33) 

min  )  d  | 

9  “ 

[P(Ef,  q)  ~ 

h). 

(34) 

i= 1 


The  functional  (32)  penalizes  absolute  differences  between  the  data  and  model.  As  illustrated  in 
Figure  9b,  it  will  produce  model  fits  which  tend  to  be  more  accurate  in  the  high  gradient  switching 
regime  than  in  the  low  gradient  saturation  region.  For  applications  which  require  high  accuracy 
throughout  the  drive  range,  two  techniques  can  be  employed  to  modify  the  manner  through  which 
the  difference  between  the  model  and  data  are  penalized.  One  alternative  is  to  employ  the  functional 
(33)  which  weights  the  data  in  the  saturation  region  more  heavily  through  a  cosine  filter  of  the  type 
illustrated  in  Figure  9c.  A  more  accurate,  but  computationally  more  intensive,  technique  employs 
the  Euclidean  metric  d  which  minizes  the  distance  between  the  model  and  data  as  illustrated  in 
Figure  9b.  For  each  functional,  initial  parameter  choices  go  are  obtained  through  the  previously 
discussed  asymptotic  relations  or  a  priori  material  information. 

4.2  Experimental  Validation  for  PZT5A,  PZT5H  and  PZT4 

To  illustrate  the  accuracy  and  flexibility  of  the  model  and  parameter  estimation  techniques  for 
a  variety  of  compounds,  we  consider  the  characterization  of  PZT5A,  PZT5H  and  PZT4  wafers. 
In  all  cases,  data  was  collected  at  200  mHz  to  minimize  frequency  effects.  For  consistency,  the 
discretized  model  (30),  with  P  given  by  (22)  or  (31),  was  employed  in  all  three  cases.  However, 
we  note  that  analogous  results  have  been  obtained  with  the  kernel  (9)  and  the  discretized  version 
of  the  statistical  mechanics  model  (29).  Finally,  the  limits  Nt  =  Nj  =  80  in  (30)  -  obtained  using 
Nq  =  Np  =  20  subdivisions  with  the  4  point  Gaussian  rule  -  ensured  convergence  of  the  Gaussian 
quadrature  routines. 


20 


PZT5A 


We  consider  first  the  characterization  of  hysteresis  exhibited  by  PZT5A  for  various  input  field 
levels.  Data  was  collected  from  a  rectangular  1.7  cm  x  0.635  cm  x  0.0381  cm  wafer  at  peak  voltages 
ranging  from  600  V  to  1600  V.  Corresponding  field  values  were  computed  using  the  relation 

E  =  V/h 

where  h  =  3.81  x  10~4  nr  denotes  the  thickness  of  the  wafer.  The  resulting  hysteretic  E-P  relations 
are  plotted  in  Figure  10. 

The  polarization  was  modeled  using  the  relation  (30)  with  the  piecewise  linear  kernel  (22).  The 
coercive  field  for  the  1600  V  data  is  1.2  x  106  V/nr  and  the  slope  after  field  reversal  in  saturation 
is  3.6  x  10s.  These  two  values  were  respectively  employed  as  initial  estimates  for  Ec  and  r/.  The 
functional  (32)  was  then  employed  to  estimate  the  parameters  Pr  =  0.04  C/nr2,  Ec  =  0.96507  x 
106  V/nr,  r)  =  9x  10s,  c  =  0.3582  V2/nr2,  b  =  2.1407  x  1011  V2/nr2  and  C  =  8.57  x  10~12  through  a  fit 
to  the  1600  V  data.  The  model  with  these  parameter  values  was  then  used  to  predict  the  hysteretic 
E-P  relation  at  the  600  V,  800  V  and  1000  V  input  levels  yielding  the  fits  plotted  in  Figure  10.  It  is 
observed  that  the  model  accurately  quantifies  the  hysteresis  through  the  switching  region  at  1600  V 
but  is  less  accurate  near  saturation  since  errors  in  this  region  are  penalized  less  by  the  absolute 
functional  (32).  The  accuracy  of  the  predictions  at  the  lower  drive  levels  attests  to  the  flexibility  of 
the  model  for  quantifying  the  E-P  relation  through  a  wide  range  of  field  inputs. 

The  use  of  the  cosine-weighted  functional  (33)  yields  the  parameters  Pr  =  0.04  C/nr2,  Ec  = 
0.866010x  106  V/nr,  r/  =  9.5xl08,  c  =  0.4272  V2/m2,  b  =  1.9754xlOn  V2/nr2  and  C  =  7.9926xl0"12 
and  produces  a  model  fit  with  slightly  improved  accuracy  in  the  saturation  region  but  less  accuracy 
near  the  coercive  field  and  at  low  drive  levels  (see  Figure  11).  Both  algorithms  provide  sufficient 
accuracy  for  quantifying  hysteresis  inherent  to  PZT5A  for  a  broad  range  of  applications. 

600  V  800  v 


Figure  10.  PZT5A  data  ( - )  and  model  predictions  ( - )  with  the  parameters  Pr  =  0.04  C/nr2, 

Ec  =  0.96507  x  106  V/nr,  rj  =  9  x  108,  c  =  0.3582  V2/m2,  b  =  2.1407  x  1011  V2/m2,  C  =  8.57  x  10"12 
determined  by  the  absolute  least  squares  functional  (32).  Abscissas:  electric  field  (MV/nr),  ordinates: 
polarization  (C/nr2). 
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600  V  800  V 


Figure  11.  PZT5A  data  ( - )  and  model  predictions  ( - )  with  the  parameters  Pr  =  0.04  C/m2, 

Ec  =  0.866010  x  106  V/rn,  rj  =  9.5  x  10s,  c  =  0.4272  V2/m2,  b  =  1.9754  x  1011  V2/m2,  C  = 
7.9926  x  10”12  determined  by  the  cosine-weighted  least  squares  functional  (33).  Abscissas:  electric 
field  (MV/rn),  ordinates:  polarization  (C/m2). 

PZT5H 

To  illustrate  the  performance  of  the  model  for  characterizing  a  second  soft  PZT  compound,  we 
consider  data  collected  from  a  3.81  cm  x  0.635  cm  x  0.0381  cm  PZT5H  wafer  at  input  levels  ranging 
from  600  V  to  2200  V.  As  with  the  PZT5A  sample,  data  was  collected  at  200  mHz  to  minimize 
frequency  effects. 

For  this  data  set,  parameters  were  estimated  through  a  fit  to  the  2200  V  data  and  the  resulting 
model  was  used  to  predict  the  hysteretic  E-P  relation  at  lower  drive  levels.  Initial  values  for  Ec  and 
r i  were  obtained  from  the  coercive  field  Ec  =  0.9  x  106  V/m  and  reciprocal  slope  4j p  =  3.7  x  108,  after 
saturation,  for  the  2200  V  data.  The  absolute  least  squares  functional  (32)  yielded  the  parameters 
PR  =  0.04  C/m2,  Ec  =  0.747690  x  106  V/m,  rj  =  6.5  x  10s,  c  =  0.2612  V2/m2,  b  =  2.8425  x 
1011  V2/m2,  C  =  1.1526  x  10-11  and  fits  depicted  in  Figure  12.  The  Euclidean  metric  (34)  yielded 
the  parameter  values  Pr  =  0.04  C/m2,  Ec  =  0.698990  x  106  V/m,  7j  =  6.5  x  108,  c  =  0.3439 
V2/m2,  b  =  3.2407  x  1011  V2/m2,  C  =  8.0932  x  10~12  and  fits  depicted  in  Figure  13.  It  is  observed 
that  because  the  absolute  metric  heavily  penalizes  discrepancies  in  high  gradient  regions,  the  fits  in 
Figure  12  are  very  accurate  in  the  coercive  region  but  saturate  too  quickly.  The  use  of  the  Euclidean 
metric  provides  more  uniform  penalties  throughout  the  drive  range  and  hence  provides  a  model  which 
is  accurate  both  in  switching  and  saturation. 

PZT4 

The  final  compound  that  we  consider  is  the  hard  material  PZT4.  Data  collected  from  a  3.81  cm 
x  0.635  cm  x  0.381  cm  wafer  at  input  levels  of  600  V  through  1800  V  is  plotted  in  Figure  14.  For 
the  1800  V  input,  the  coercive  field  1.42  x  106  V/m  reflects  that  more  energy  is  required  to  turn 
dipoles  than  in  the  softer  PZT5  compounds. 
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Figure  12.  PZT5H  data  ( - )  and  model  predictions  ( - )  with  the  absolute  least  squares 

functional  (32)  used  to  determine  the  parameters  Pr  =  0.04  C/m2,  Ec  =  0.747690  x  106  V/m, 
rj  =  6.5  x  10s,  c  =  0.2612  V2/m2,  b  =  2.8425  x  1011  V2/m2,  C  =  1.1526  x  10“12.  Abscissas:  electric 
field  (MV/rn),  ordinates:  polarization  (C/m2). 


Figure  13.  PZT5H  data  ( - )  and  model  predictions  ( - )  with  the  Euclidean  least  squares 

functional  (34)  used  to  determine  the  parameters  Pr  =  0.04  C/m2,  Ec  =  0.698990  x  106  V/m, 
t /  =  6.5  x  10s,  c  =  0.3439  V2/m2,  b  =  3.2407  x  1011  V2/m2,  C  =  8.0932  x  10~12.  Abscissas:  electric 
field  (MV/rn),  ordinates:  polarization  (C/m2). 
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600  V  1000  V 


Figure  14.  PZT4  data  ( - )  and  model  predictions  ( - )  with  the  parameters  Pr  =  0.045  C/m2, 

Ec  =  1.05  x  105  V/rn,  i]  =  4.0  x  108,  c  =  0.3018  V2/m2,  b  =  2.1924  x  1011  V2/m2,  C  =  6.8287  x  10~12. 
Abscissas:  electric  field  (MV/m),  ordinates:  polarization  (C/m2). 


The  three  functionals  (32)-(34)  yield  roughly  equivalent  parameters  and  the  model  response  with 
the  parameters  Pr  =  0.045  C/m2,  Ec  =  1.05  x  105  V/m,  r]  =  4.0  x  108,  c  =  0.3018  V2/m2, 
b  =  2.1924  x  1011  V2/m2,  C  =  6.8287  x  10~12,  estimated  through  a  fit  to  the  1800  V  data,  is 
compared  with  the  data  in  Figure  14.  It  is  observed  that  while  model  is  very  accurate  in  the  high 
drive  regime,  it  does  not  fully  quantify  the  sharp  transition  prior  to  coercivity  due  to  limitations 
resulting  from  the  choice  of  the  lognormal  and  normal  functions  used  to  respectively  characterize 
the  densities  of  the  coercive  and  effective  fields.  This  tendency  is  accentuated  when  the  model  is 
subsequently  used  to  predict  the  600  V  E-P  relation.  While  the  model  provides  sufficient  accuracy 
for  applications  such  as  control  design,  the  discrepancy  illustrates  that  improvements  can  be  made 
in  the  model  when  modeling  hard  compounds.  We  are  presently  investigating  the  development  of 
higher-order  kernels  and  techniques  for  identifying  general  densities  to  provide  additional  accuracy 
for  high  performance  applications  which  utilize  PZT4  transducers. 

4.3  Biased  Asymmetric  Minor  loops 

To  illustrate  the  capability  of  the  model  to  quantify  biased  minor  loops,  the  field  plotted  in 
Figure  15a  was  input  to  the  model  to  generate  the  E-P  relation  depicted  in  Figure  15b.  The 
parameters  were  taken  to  be  Pr  =  0.04  C/m2,  Ec  =  0.698990  x  106  V/m,  r/  =  6.5  x  108,  c  = 
0.6877  V2/nr2,  b  =  3.2407  x  1011  V2/m2,  C  =  8.0932  x  10~12  which  are  close  to  those  identified 
for  PZT5H  (the  parameters  were  modified  slightly  to  highlight  aspects  of  the  biased  loops).  Loop  1 
illustrates  the  capability  of  the  model  to  characterize  the  quadratic  Raleigh  behavior  experimentally 
observed  for  low  drive  levels  whereas  Loop  2  illustrates  that  the  model  enforces  closure  of  biased, 
asymmetric  minor  loops  on  the  initial  polarization  curve.  The  continuity  in  slope  of  the  initial 
curve,  following  the  closure  of  Loop  2,  illustrates  that  the  model  incorporates  the  deletion  property 
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(a)  (b) 

Figure  15.  (a)  Field  input  E,  and  (b)  polarization  predicted  by  the  discretized  model  (30). 


which,  along  with  congruency,  forms  one  of  the  necessary  and  sufficient  conditions  for  classical 
Preisach  models  [8] .  Loops  3  and  4  illustrate  the  ability  of  the  model  to  enforce  closure  of  multiply 
nested  loops  while  Loop  5  encapsulates  biased  behavior  near  saturation.  When  combined  with  the 
experimental  results  in  Section  4.2,  the  behavior  depicted  here  provides  the  model  with  substantial 
flexibility  for  material  characterization  and  control  design. 

5  Concluding  Remarks 

The  theory  presented  in  this  paper  provides  a  technique  for  quantifying  constitutive  nonlinearities 
and  hysteresis  inherent  to  piezoelectric  compounds  in  a  manner  conducive  to  bulk  material  char¬ 
acterization  and  model-based  control  design.  Through  the  combination  of  free  energy  analysis  at 
the  lattice  or  domain  level  with  the  assumption  that  certain  physical  parameters  are  stochastically 
distributed,  macroscopic  models  having  a  small  number  of  parameters  (5-6)  are  constructed.  Fur¬ 
thermore,  several  of  these  parameters  can  be  correlated  with  physical  attributes  of  the  E-P  data 
to  facilitate  parameter  identification  and  model  updating  to  accommodate  changing  operating  con¬ 
ditions.  The  model  accommodates  transient  dynamics  in  the  E-P  relation,  enforces  return  point 
memory,  and  guarantees  the  closure  of  both  symmetric  and  biased,  asymmetric  minor  loops.  The 
model  does  not  enforce  congruency  in  the  saturation  regions  of  the  E-P  curve  which  reflects  the 
measured  behavior  of  the  materials  in  these  regions. 

The  numerical  examples  and  comparison  with  experimental  data  illustrate  low-frequency,  fixed- 
temperature  capabilities  of  the  model.  In  its  present  formulation,  however,  it  also  incorporates  certain 
relaxation  mechanisms  and  temperature-dependencies.  The  latter  effect  has  been  validated  in  the 
context  of  relaxor  ferroelectric  compounds  (PMN)  and  a  comprehensive  set  of  examples  illustrating 
the  quantification  of  these  properties  in  PZT  will  appear  in  a  future  publication.  The  model  in  its 
present  form  does  not  incorporate  polarization  changes  due  to  variable  stresses  and  these  extensions 
are  under  current  investigation. 

The  use  of  Algorithm  1  renders  the  model  highly  efficient  so  that  multiloop  simulations  run  in 
the  order  of  seconds  on  an  aging  and  overworked  workstation.  In  the  algorithm  described  in  [34],  the 
nronotonicity  of  the  modeled  E-P  relation  was  invoked  to  construct  an  inverse  model  for  inclusion  in 
feedforward  or  feedback  loops  to  linearize  (at  least  approximately)  the  response  of  actuators  operating 
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in  nonlinear  and  hysteretic  and  nonlinear  regimes.  The  analogous  magnetic  model  and  inverse  are 
employed  in  [17]  to  construct  robust  control  designs  to  achieve  high  accuracy,  high  speed  position 
control  for  Terfenol-D  transducers.  Hence  the  unified  nature  of  the  modeling  approach  facilitates 
unified  control  designs  that  utilize  model-based  inverse  compensators. 
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